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Abstract 
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blade chord 


This paper discusses the numerical simulation M 
of three-dimensional transonic flow about propel- 
ler blades. The equation for unsteady potential 
flow about propellers is given for an arbitrary Mr 

coordinate system. The small disturbance form of 
the equation is then derived for a new helical 
coordinate system. The new coordinate system is Mj 

suited to propeller flow and allows cascade bound- 
ary conditions to be straightforwardly applied. 

A numerical scheme is employed which solves the N 

steady flow as an asymptotic limit of unsteady 
flow. Solutions are presented for subsonic and Nr 

transonic flow about a 5 percent thick bicircular 
arc blade of an eight bladed cascade. Both high q' 

and low advance ratio cases are given, and include 
a lifting case as well as nonlifting cases. The R 

nonlifting cases are compared to solutions from 
a Euler code. r 


Nomenclature 

t 

A function relating periodic helical 

coordinates to Cartesian t 

coordinates 

tr 

A 1 * A 2 * - • • * a 8 coefficients of small disturbance 

equation, Eq. (51) U 


a undisturbed fluid velocity vector u',v',w‘ 

in general rotating coordinate 
system _ 

V 

a 0 nonnegative constant, Eq. (38) 

B function relating k to z, 

Eq. (32) X , Y, Z 


C 00 

n Y 

E 

F 


speed of sound at infinity 

x,y,z 

special difference operator in 
mean flow direction 

xl,x2, x 3 

damping function in helical 
coordinate system 

y 1 ,y z ,y 3 

function defining blade surface 


F 0 ,F U lower and upper blade functions Greek : 


g 

9ij,g ij 

h 

ho»h u 

k 


determinant of the metric tensor 

covariant and contravari ant com- 
ponents of the metric tensor 

blade profile parameter 

lower and upper blade profiles 

ratio of specific heats 


a 

T 

Y,r,s 

Vr,T 


a 

6 


^Member AIAA. 


radially dependent Mach number 
of free-stream 

Mach number of propeller tip at 
r = R 

circumferential Mach number of 
propeller tip at r « R 

number of ADI iterations 

number of blades in cascade 

perturbation velocity 

radius of propeller 

radius from axis of rotation; 
second helical coordinate 

radius from origin 

time 

maximum thickness of blade 
helical velocity 

perturbation fluid velocities in 
x,y,z-directions 

advance velocity of propeller; 
relative velocity between iner- 
tial and rotating frames 

Cartesian coordinates in inertial 
frame 

Cartesian coordinates in rotating 
frame 

Cartesian coordinates correspond- 
ing to x,y,z 

helical coordinates corresponding 
to y,r,e 


angle of attack 

circulation function 

helical coordinates 

uniform computational coordinates 

increment operator 

ratio of maximum blade thickness 
to chord length, 6 = t max /e 



6 r t; finite difference operators 

e inverse of blade aspect ratio, 

€ = i/R 

e-j streamwise switching operator 

permutation symbol 

0 circumferential angle about axis 

x advance ratio, x = V/(rc R) 

C m half the distance between neigh- 

boring blades 

T s as measured from nearest 

blade 

$ total velocity potential 

cp perturbation velocity potential 

n rotational speed of propeller 

a) frequency of airfoil motion 

Subscripts and superscripts : 

I , J,K nodal indices of computational 

mesh 

i,j,k indices for tensor quantities 

L.E. leading edge 

n,n+l value at time step n and n+1 

T.E. trailing edge 

* value at first intermediate time 

step 

** value at second intermediate 

time step 

1 derivative of quantity 


Introduction 

Computational fluid dynamics (CFD) has become 
a powerful tool in the analysis of complex fluid 
flows. Coinciding with the arrival of increasingly 
powerful CFD methods, renewed interest has emerged 
in propellers which offer large performance 
improvements for aircraft that cruise in the high 
subsonic speed range. This paper presents results 
of the numerical simulation of the flow about pro- 
peller blades operating with tip speeds in the 
transonic regime. These results were obtained 
using an ADI finite difference procedure to solve 
the low frequency form of the three-dimensional 
small disturbance equation. Expressing the small 
disturbance equation in terms of novel helical 
coordinates significantly facilitates the numeri- 
cal solution for propeller flow. 

Recent propeller research programs have 
introduced turboprops which differ considerably 
from their low-speed turboprop counterparts. The 
most notable changes are the use of eight to ten 
thin, highly swept, small -diameter blades instead 


of two to four longer, thicker and unswept conven- 
tional blades. These changes allow propellers to 
operate efficiently in the transonic flow regime. 
However, they require sophisticated computational 
methods in order to predict and calculate their 
aerodynamic characteristics. This is mainly due 
to the increased three dimensionality of the flow 
field encompassing these low aspect ratio blades 
and to the rise of mathematical nonlinearity in 
the flow model which is associated with their 
transonic tip speeds. Also, their thin structural 
design makes them prone to flutter. This is a 
source of unsteadiness which must be considered in 
addition to the inherent unsteadiness of the tran- 
sonic regime. 

Many examples of both steady and unsteady 
flows of practical interest are well approximated 
by various forms of the governing equations of 
fluid mechanics. The recent marked advance in 
computer capabilities and the rapid growth of 
efficient and accurate numerical methods have 
allowed a number of previously untenable problems 
to be examined. The use of numerical methods 
offers the the ability to handle complex equations 
with intricate geometries such as exist in the 
case of simulating flow about rotating propeller 
blades. However, as with any method, limits exist 
on the capabilities. These limits involve the 
speed and memory size of the computers and the 
availability of efficient algorithms as well as 
considerations of robustness and affordability. 

For a given problem, judgement dictates the trade- 
off between the suitability of a particular model 
and the feasibility of the calculation. A number 
of numerical codes are currently being developed 
to calculate propeller flow. These solve either 
the full potential or Euler equations. Although 
these equations accurately model the flow, they 
are computationally more expensive and more com- 
plicated to numerically code than the small dis- 
turbance equation. 

A numerical code was developed*- 4 to simulate 
the flow over blade-tips of helicopter rotors. 

This code solves the small disturbance equation 5 
appropriate to a helicopter in forward flight. In 
particular, results were obtained that showed the 
code was able to track the development of a shock 
and its subsequent propagation upstream. The suc- 
cessful application of the small disturbance code 
in this manner indicated that the same approach 
could be applied to the propeller problem. 

The fact that turboprops operate at high sub- 
sonic cruise speeds means that at least the tip 
region of each blade will be embedded in a tran- 
sonic helical flow. Shocks may likely exist ahead 
of or on the blades, depending on the flight speed 
and propeller design. Shocks of high strength 
indicate large shock wave drag. Efficient flight 
dictates that only weak shocks are acceptable, 
which implies little entropy generation. The 
potential formulation should provide a good approx- 
imation to the Euler equations as long as the 
shock strengths are less than that of a normal 
shock with an upstream Mach number less than about 
1.3. Presently, turboprops that have helical Mach 
numbers slightly less than one at the hub to about 
1.12 at the tip are being designed. It is reason- 
able to expect that the shock strength is suffi- 
ciently weak near the blade tips, where three- 
dimensional relief effects occur, so that the 
potential formulation is adequate. Furthermore, 
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the use of a perturbation potential in place of a 
full potential significantly simplifies the formu- 
lation of the problem. 

Although many similarities exist between the 
equations developed for the helicopter rotor and 
the propeller problem, fundamental differences 
exist between the two. These differences demand 
special consideration be made in accurately treat- 
ing each type of flow. For example, in the case of 
the helicopter rotor, the flow field as observed 
in a coordinate system attached to the blade 
remains unsteady even if blade flutter is ignored. 
However, when the axis of rotation for a propeller 
coincides with its flight direction, a steady flow 
results with a blade-attached coordinate system. 
This allows the steady problem to be examined 
separately from the flutter or gust problem. With 
the goal of using the numerical algorithm devel- 
oped for the rotor code to the propeller problem, 
several concerns need to be addressed. Many of 
these concerns are naturally handled by using a 
suitable coordinate system. 

The key to simplifying the potential formula- 
tion of the propeller problem is the use of the 
helical free-stream direction as the primary direc- 
tion. This provides an accurate primary flow upon 
which a perturbation can be superposed. The use 
of a helical flow, which inherently captures the 
flow curvature, results in a marked improvement 
over using an axial flow as the primary flow. The 
optimum flow to base a perturbation about is the 
exact flow; helical flow much more closely approxi- 
mates the exact flow about a rotating and advanc- 
ing propeller than does axial flow. 

Since the helical flow direction captures the 
fundamental properties of propeller flow, rela- 
tively few terms are needed in the resulting dis- 
turbance equation to provide accurate modeling of 
the flow. In fact, no more terms are needed than 
when using Cartesian coordinates for the rotor 
problem. In addition, due to the similar nature 
of the two problems, the terms in the two sets of 
disturbance equations correspond identically except 
for their coefficients. Thus, the propeller equa- 
tion is amenable to the same potential solver as 
developed for the helicopter code. This is quite 
advantageous since much effort has been expended, 
both in making the potential solver efficient and 
in verifying its operation. 

While the helical flow direction is an ideal 
choice for the primary flow direction, other con- 
siderations must be taken into account in order 
to provide a satisfactory coordinate system. These 
considerations hinge on the desire to include cas- 
cade effects and the employment of small disturb- 
ance boundary conditions near the blade surfaces. 
Proper specification of the helical coordinates 
permits the straightforward treatment of annular 
cascade flows while maintaining orthogonal proper- 
ties near the blade locations. The helical coor- 
dinates given in this paper allow both of these 
requirements to be met. 

These helical coordinates are analytically 
defined in terms of basic propeller parameters 
through simple transformations to Cartesian coor- 
dinates. Thus, they automatically adapt to the 
cascade configuration for any number of blades 
and blade twist, which itself is a function of 


the propeller advance ratio. The helical coordi- 
nates enter the problem in terms of the metric 
tensor. The components of the metric tensor are 
readily evaluated in a separate subroutine. The 
Jacobian of the transformation is particularly 
simple and provides a means to determine if a 
given transformation is valid. 

Obtaining accurate solutions for steady, yet 
alone unsteady, three-dimensional transonic flows 
presents a difficult challenge even for the most 
advanced computational methods. The present work 
assumes a steady state flow exists as observed in 
a reference frame attached to the propeller system 
which is rotating with a constant angular velocity. 
This is considered a necessary step to be com- 
pleted before undertaking the unsteady problem. 

In addition, the fluid will be regarded as being 
inviscid in an irrotational flow field which is 
expressible in terms of a potential formulation. 

The solution is sought in terms of a disturbance 
potential. This disturbance potential is the 
potential associated with the reduced velocity 
obtained by subtracting both the free-stream 
velocity and the rotational velocity, resulting 
from the transformation to a noninertial reference 
frame, from the total velocity. This solution 
should be valid for lightly loaded blades. 

This paper presents the governing equation 
for the three-dimensional small perturbation 
potential and its derivation from the potential 
equation valid for a rotating reference frame. 

This is given in general tensor form applicable 
to any suitable coordinate system, which need not 
be orthogonal. An original helical coordinate 
system is then specified which aligns with the 
mean undisturbed flow and provides for easy 
application of periodic boundary conditions as 
well as providing appropriate surfaces on which 
to enforce the small disturbance boundary 
conditions. The small disturbance form of the 
equation is then presented in terms of these 
coordinates which are scaled by the relevant 
transonic scaling parameters. A generalized form 
of the Douglas-Gunn algorithm is given which 
allows for additional cross-derivati ve terms to 
be included in the potential solver. However, 
for the flow cases presented here, these terms 
had little effect on the solution, and, hence, 
for the results of this effort they were omitted 
from the governing equation. 

Solutions are presented for an annular 
cascade consisting of eight bicircular arc blades 
with a maximum thickness of 5 percent. These 
blades are run for two geometric operating 
conditions: a high advance ratio case where the 
advancing speed is large compared to the rotational 
speed of the propeller tip, and a low advance 
ratio case where these two speeds are equal. For 
the high advance ratio case, the helical free-stream 
Mach number relative to the blade tip is Mr =* O.B. 
For the low advance ratio case, results for two 
values of the free-stream Mach number are 
presented - a subsonic free-stream case of 
Mr = 0.8, and a transonic case of Mr = 1.1 
which is representative of a turboprop. In 
general, the results are for zero angle of attack; 
however, a lifting case is presented for the low 
advance ratio and high free-stream condition. The 
nonlifting cases are compared with results from an 
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Euler code. In the discussion and results sections, 
the abbreviation HSD will be used to denote the 
helical small disturbance computations. 

Governing Equation 

This section presents the potential equation 
and the development of an approximation to it 
which is valid for small disturbance flow about 
thin propeller blades which are lightly loaded. 
First, the potential equation in rotating coor- 
dinates will be given as developed for a coordinate 
system attached to a helicopter blade. Following 
this, the equivalent tensor equation valid for any 
curvilinear coordinate system will be presented 
for an accelerated system. Finally, the small 
perturbation equation for generalized helical 
coordinates will be given. 

Potential Equation in Noninertial System 

The potential equation can be expressed in 
vector form° as 

3 2 1 2 

+ -^^g) * 74*7(7$) s= 

cf - (t - 1) |* t + |(v») 2 j (1) 

where $ is the velocity potential, c ro is the 
upstream speed of sound and 4 is the ratio of 
specific heats. 

The potential equation can be expressed in 
an xyz coordinate system which is both trans- 
lating with an arbitrary constant velocity V 
away from the inertial frame, and rotating with_ 
an arbitrary constant rotational velocity 0 x i . 
The angular velocity vector n (with magnitude 
n) has the direction of the axis through the_ 
origin about which xyz are rotating, and 7 is 
the position vector in this system. The relation- 
ship between the two coordinate systems is shown 
in Fig. 1. Although both V and IT are general 
vectors, they will soon be restricted to a common 
axial direction. 

A potential equation can also be established 
in the noninertial frame. This is initiated by 
defining a perturbation potential <p which 
separates the contribution of the free-stream from 
the total potential. The free-stream contains 
both the uniform flow and the bulk rotational 
flow, and 


<P tt + a*v(v<p*a) + 2a*V(o t - pxa*vcp 


3 V — 2 

- -^jr*7<p + 2vcp*v<p t + a*v(v<p) 

= j c f “ (^“1) ^ + a*7cp + -j-(vcr) jv^<p (4) 

where a is the negative of the vector sum of the 

translational velocity and the rotational velocity 
as seen by an observer_in the rotating frame. Thus, 
for such an observer, a can be simply written as 

a = - (V + H. x 7) (5) 

It should be pointed out that in the trans- 
formation from the potential Eq. (1) for an iner- 
tial reference frame to Eq. (4) for a noninertial 
frame, the cubic term, the last term on the left- 
hand side of Eq. (1), was dropped. Otherwise, 

Eq. (4) remains an exact equation. The introduc- 
tion of the perturbation potential does not, in 
itself, imply any approximations. 

Potential Equation in Noninertial Curvilinear 
Coordinate System 

The application of Eq. (4) to the propeller 
problem can be carried out by further simplifica- 
tion in much the same manner as was done in deriv- 
ing the small disturbance equations for the 
helicopter problem. However, since the helical 
coordinate systems which will be used for the pro- 
peller problem are nonorthogonal , the potential 
Eq. (4) is more tractable if expressed in invariant 
tensor form. The tensor form of Eq. (4) is 



q* = Vq> = q - (V + ox*) (2) 

represents the irrotational perturbation velocity. 
The gradient operator is invariant with respect to 
choice of Cartesian coordinate systems. The com- 
ponents of the perturbation velocity 'q' are 
defined as 


u ' = 


o cp 

3 X ’ 


3(P 

ay* 


3 cp 
3 Z 


(3) 


The transformation of the potential equation has 
been carried out by Isonr for a coordinate system 
translating and rotating with constant, but 
otherwise arbitrary, V and p. The potential 
equation in the noninertial system is given as 


for a translational velocity V which is independ- 
ent of time. In this equation ai and aj 
represent the contravari ant and covariant compon- 
ents of the general velocity .vector 7 defined by 
Eq. (5). The quantities g^ are the contra- 
variant components of the metric tensor for the 
transformation between the curvilinear coordinates 
yi and our orthogonal Cartesian coordinates. The 
quantity g is the determinant of the correspond- 
ing covariant components of the metric tensor g-j -j 
which do not appear explicitly in the equation. 

The symbol is the permutation symbol which 

equals zero for repeating values of i,j,k, unity 
for cyclic (even) permutations of 1,2,3, and neg- 
ative unity otherwise. 
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Potential Equation in Helical Coordinates 


Equation (6), expressed in terms of a rotating 
coordinate system, can be simplified considerably 
by choosing one of the curvilinear coordinates to 
be in the direction of the vector a. From Eq. (5), 
this is seen to be the direction opposite the 
vector sum of the translational velocity and rota- 
tional velocity which we now take to be orthogonal 
such that V is aligned with the axis of rotation, 
as shown in Fig. 2. With this arrangement the 
vector a is in the direction of the free-stream 
velocity vector as it appears to an observer whose 
frame of reference rotates with the blade. In this 
instance the magnitude of IT will be designated 
by U, such that in matrix form 


a 


0 


(7) 


applied to Eq. (6) for the case of a propeller 
blade using helical coordinates. The use of hel- 
ical coordinates introduces the metric tensor into 
the approximate equation and allows representation 
of flow curvature. 

The derivation of the approximate equations 
for flow about a propeller is based on considera- 
tion of the following parameters: the free-stream 

axial Mach number M*,; the thickness ratio 6; the 
ratio of the chord i to the blade-tip radius R, 
which is the inverse of the aspect ratio, e; and 
the advance ratio x. For true unsteady problems 
the reduced frequency would enter; however, since 
we are examining only steady solutions, the reduced 
frequency need not be introduced. Also, for cas- 
cade solutions, a parameter representing the block- 
age of the flow, say solidity, would normally enter 
here the blades will be considered far enough apart 
that variations in solidity will not be important. 
For an advanced turboprop, typical values of the 
parameters are = 0.8, 6 = 0.02, e « 0.3, and 
x = 1. 


and 


U 2 = (nr) 2 + V 2 (8) 

where r is the radial distance as measured from 
the axis of rotation. Equation (6) then simpli- 
fies to 



The lateral direction y 3 , which lies essen- 
tially normal to the blade, is scaled by dividing 
it by the value a 3 ' 3 . This could be done for the 
radial-like or spanwise direction y 2 , also since 
the physical justification for scaling is to trans- 
form the lateral dimensions to account for the 
weak diminution of disturbances in these directions 
in comparison to the streamwise direction. How- 
ever, Isom chose to scale the spanwise direction 
differently by introducing e. This is accom- 
plished by normalizing the steamwise direction 
yl along with the already scaled y 3 direction 
by while normalizing y 2 .by R. Since these 
two methods of scaling the y ^ result in the 
same final approximation, either can be used. The 
latter choice is made here and the directions are 
scaled such that the original coordinates trans- 
form to the dimensionless ones, as given by 


L 3y3r yg 3 y \ / arj 


(9) 


Approximate Potential Equation in Scaled Helical 
Coordinates 


Equation (9) is exact in that it is equivalent 
to Eq. (4), but expressed in a blade-fixed refer- 
ence system. However, this equation is much too 
complicated to be solved efficiently. Therefore, 
a systematic simplification is necessary to arrive 
at an approximation which retains the nonlinear 
features of the problem and which is valid for 
small disturbances about the mean flow. Starting 
from Eq. (4), such a procedure was carried out by 
lson£ using Cartesian coordinates for flow about 
the tip region of a helicopter blade. Following 
the development of the approximate equation for 
the helicopter problem, a similar process is 


y 1 = ^y 1 , y 2 = Ry 2 , y 3 = -jjj y 3 (io) 

where the tilde denotes the dimensionless coordi- 
nates. In addition, time and the disturbance 
potential are nondimensional ized as follows: 

t = ^ , to = <p (11) 

The use of the tilde is for clarif ication only, 
and in the following equations it will not be 
used. 

By defining 

M = (12) 
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and introducing a constant Mach number character- 
istic of the rotational speed of the propeller 
tip which is defined by 

(13) 

oo 

the perturbation equation can be derived and is 
given in scaled coordinates as 


where 

F = F u (y 1 ,y 2 ,y 3 ,t) = y 3 - hJyV.t) = 0 (18a) 

on the upper (suction) blade surface from the lead 
ing edge (L.E.) to the trailing edge (T.E.), and 

F = F q ( y 1 , y 2 , y 3 , t ) = y 3 - h^.y 2 ,!) = 0 (18b) 
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For low frequency or steady problems, it is 
admissible to delete the second derivative with 
respect to time. This gives the final form of 
the equation in terms of the scaled variables as 
being 
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(15) 


on the lower (pressure) blade surface. Here, h is 
the profile parameter. This follows the convention 
used for wing surfaces. The coordinates are aligned 
as shown in Fig. 3, where y 1 is nearly along the 
mean chord, y 3 is nearly normal to the chord, and 
y 2 lies in the spanwise direction. 

One of the key developments of this work is 
the generation of coordinate systems possessing 
the characteristic that g 13 ~ 0 at the blade 
surfaces. For these coordinate systems, the 
unsteady boundary condition at the blade surface 
is simplified. 

The general time-dependent boundary condition 
which allows the profile parameter to include 
pitching, bending, and twisting is given in scaled 
variables as 




dh 

at 


M ah 
M. 


T ay 


(19) 


This time dependence is included mainly for any 
future efforts which consider flutter. It is 
interesting to note that in dimensional Cartesian 
coordinates (y 1 = x, y 2 = y, y 3 = z), this is the 
f amil iar 


a? ^ ah ^ y a h 
7z ~ 1ft u ax 


( 20 ) 


Finally, the steady scaled boundary condition 
is simply 


3a, 1 M ah 

77 7 s ' "T 77 


(a) 


Equation (15) is the small disturbance equation 
which is solved numerically by the ADI Douglas- 
Gunn algorithm. 

Boundary Conditions 


This boundary condition is not actually applied 
at the airfoil surfaces but rather at two surfaces 
separated slightly such that they straddle the 
mean chord surface which itself lies on, or near, 
the y 3 = 0 surface. Thus, at the upper surface 


Wall Tanqency Condition for Blade Surfaces 

The small disturbance boundary condition to be 
applied at the blade surfaces will now be estab- 
lished. This necessitates returning to the 
unsealed variables momentarily. 


3a> 



y 3 =0+ 


1 M 
.33 M 


3h, 
T ay 


and, at the lower surface 


(22a) 


In blade-fixed coordinates, if the equation 
for the surface of a blade moving in a time- 
dependent manner is 

F(y 1 ,y 2 ,y 3 ,t) = o, (16) 


3(p 

17 


y 3 =0- 


1 M 
7™ 


alv 


T ay 


Flow Field Boundary Conditions 


(22b) 


then the vanishing of the fluid velocity component The U p S tream boundary condition on the pertur- 

normal to the surface brings bation potential is 


F ( y 1 , y 2 , y 3 , t ) = ^ + (a + v®)-vF = 0 (17) 


(p = 0 


(23) 
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This merely states that the free-stream condition 
prevails sufficiently far upstream. 

In the downstream direction, the perturbation 
potential is not known, but at large distances the 
condition 



(24) 


Periodic Helical Coordinate System for Cascades 

This section contains the development of spe- 
cial periodic helical coordinates suitable for 
propeller problems which include cascade effects. 
The coordinate transformation between this system 
and an orthogonal Cartesian system is specified 
so as to provide simple periodic boundaries, and 
also orthogonal properties at the blades surfaces. 


is reasonable and is used as the downstream bound- 
ary condition. 

At the hub or inner spanwise boundary, the 
condition 



(25) 


is used, and the outer spanwise boundary condition 
is simply 


<p = 0 (26) 

Although y^ is not strictly a radial direction, 
it closely approximates it for the coordinate 
systems employed here. 


Above and below the blade, cascade-type peri- 
odic boundary conditions are applied such that 


<p(y 1 ,y 2 .t) 


3_ upper 
* ^surface 


is satisfied over the boundary surfaces which are 
chosen to lie at the midchannel between the blades. 
This, in return for straightforward application, 
requires a coordinate system which provides nodal 
coordinates having azimuthal periodicity. Such 
periodic coordinate systems are discussed in 
another section. 


It remains to discuss the boundary conditions 
in the vortex wake. Rollup of a wake is ignored 
and, furthermore, the resulting vortex sheet is 
approximated as lying between the downstream 
extensions of the helical surfaces upon which the 
blade boundary conditions are imposed. The 
strength of the vortex sheet is assumed to be 
preserved as it convects downstream. In addition, 
the direction of the vortex vector is assumed to 
be "parallel 0 to the free-stream direction. Since 
the vortex sheet is a free surface, it cannot sup- 
port a pressure difference. The boundary condition 
at the trailing edge can be written in terms of 
the scaled coordinates as 


The helical coordinate system used here pro- 
vides natural surfaces in the form of helicoids 
which pass through the mean blade locations. This 
gives convenient surfaces essential to the direct 
application of blade boundary conditions. Pre- 
scribing a second coordinate so as to make it 
intersect the mean blade surfaces orthogonally, 
and yet be axially periodic in the circumferential 
direction, allows both accurate specification of 
blade boundary conditions and straightforward 
application of blade periodicity. A third coor- 
dinate lies essentially in the radial direction. 

The helical coordinate system is thereby composed 
of a coordinate aligned with the mean flow direc- 
tion, a radial-like coordinate and a coordinate 
in a circumferential-like direction. Thus, only 
in the case of the first coordinate direction are 
the coordinate curves helices. An example of the 
coordinate system in relationship to a Cartesian 
one is shown in Fig. 4 

The main advantage of these coordinates is 
that the undisturbed flow direction can be aligned 
with the helical coordinate direction y 1 . These 
streamwise helices are constructed such that the 
helices of greater radii have smaller advance 
ratios, i.e., the smaller the radius the steeper 
the spiral, exactly characterizing the mean pro- 
peller flow. 

To provide circumferential periodicity, the 
y 3 coordinate must be constructed such that, in 
a periodic fashion, it repeats its axial locations. 
The simplest choice is to make y 3 the circumfer- 
ential direction, forcing it to be independent of 
the axial value, i.e, identical to the angular 
coordinate of the familiar circular cylindrical 
coordinate system. The drawback in doing this is 
that the y 3 direction would no longer be normal 
to the streamwise direction. Orthogonality in these 
two directions is an important quality, especially 
near the surfaces of the blades. To provide 
orthogonality at the mean blade locations, the 
y 3 direction must be perpendicular to the stream- 
wise helices, at least near the blades. Away from 
the blades, the y 3 coordinate should reverse its 
axial direction so as to bend back and regain its 
original axial station. This must be repeated in 
a periodic fashion about the axis so as to conform 
to the locations of the blades. 


W T.E. = r(y2,t) (28) 

where r is the circulation around the propeller 
blade, and 

[<p] = (cp) ^ - (cp) ^ (29) 

y =0+ y =0- 

is the jump in the perturbation potential across 
the wake. 


The coordinate direction y2 which serves to 
measure the radial value is not, in general, a 
straight line. This coordinate direction is only 
a straight line at its coordinate axis, which was 
chosen to align with the leading edge of one blade, 
and at the leading edge locations of the remaining 
blades in the cascade and at lines midway between 
the leading edges. Only at the leading edge is 
the coordinate system truly orthogonal. A notable 
characteristic is that, for high advance ratios, 
they tend toward orthogonal cylindrical coordi- 
nates. This can be easily seen by inspecting the 
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metric tensors for this transformation which are 
provided in Appendix A 7 , 


z = - y l + A(r)BU) (34) 


A set of helical coordinates which has these 
properties is given by the fol lowing .transforma- 
tion to the Cartesian coordinates x 1 . 


xl = y 2 sin e 

(30) 

x 2 = y 2 cos e 

(31) 

X 3 = - -Jy 1 + A(y 2 )B(y 3 ) 

(32) 


where e is the circumferential angle measured 
as shown in Fig. 4, and which is related to the 
helical coordinates by 


B 


£ y 1 

u y 



(33) 


The total helical velocity U is, of course, a 
function of the radius. 


It is convenient to make the following assign- 
ments which are consistent with Appendix A of:' 
x 1 = x, x 2 = y, x 3 = z for the Cartesian coordi- 
nates and yl = y, y 2 = r, y^ = z for the helical 
coordinates. Using these replacements the next 
section will discuss the proper choices for A(r) 
and B(s). 

Tailoring the Periodic Helical Coordinates 

The functions A(r) and B(c) are used to 
tailor the helical coordinate system to two sets 
of the streamwise helical sheets. The two sets 
of sheets of constant z are illustrated in 
Fig. 5. They are evenly spaced in the circumfer- 
ential direction. The first set of these sheets 
will contain the mean position of uncambered, 
symmetric, twisted blades. These correspond to 
the advance helicoids. The second set of helical 
surfaces are similar and are chosen such that each 
sheet lies midway between two neighbors of the 
first set, thus forming an alternating arrangement 
of periodically spaced helicoidal surfaces. The 
boundary conditions for the airfoil surfaces are 
applied very near members of the first set, 
whereas the periodic conditions are enforced on 
the second. It should be mentioned that for asym- 
metrical blades or blades with camber, the mean 
blade position will not quite coincide with a 
helical sheet. The mean blade positions are 
assumed to lie near the first set of helical 
sheets so that small disturbance boundary condi- 
tions can be accurately applied. 


This specifies the axial coordinate in terms of 
the helical variables. It will be shown that, by 
choosing suitable forms for A and B, two of the 
metric tensors gi 3 and g ?3 become equal to 
zero concurrently at special values of z. The 
y and z curves will be orthogonal where gi 3 = 0 
and, likewise, the r and z curves will be 
orthogonal where g ?3 * 0. When these are both 
zero, the z curves will be normal to the sur- 
faces of constant z. The first objective is to 
arrange this to occur for those helical sheets 
containing the mean positions of the propeller 
blades. To assure that these surfaces are period- 
ically spaced in the circumferential direction, 

B(0 is expressed as the product of two functions 
of z in the form 

B(5) = E(0 ^ sin 4-» > \ (35) 


where E is a damping function and z m is a 
positive constant establishing the period and is 
exactly half the distance between the blades. 

Meeting the conditions of orthogonality 
requires the inspection of the functional depend- 
ence of gi 3 and g23* This will reveal what the 
final forms of A and B must be. First, the 
metric tensor g ^3 will be discussed, and then 
g 23 where additional constraints will be imposed 
on the axial velocity V(r). 

From Appendix A, 7 we have 

9 13 “ITf-II A(r)B,(F) (36) 

where the prime indicates differentiation with 
respect to the indicated argument. By using 
Eq. (35), this becomes 



+ E(r) cos irjj (37) 

As mentioned, E is a damping function and is 
needed to assure a valid coordinate transformation. 
A useful form is E = exp(-a 0 i s/s m i ) . The signi- 
ficance of the over-bar above z in the exponen- 
tial argument will be discussed later. 


As mentioned, the helical sheets are surfaces 
of constant s. It is desirable for satisfactory 
application of blade surface boundary conditions 
that the e-coordinate be orthogonal to the mean 
blade surface. In addition, simple handling of 
periodicity requires that there be no net change in 
axial distance when traversing a e-coordinate 
line from one periodic boundary to another. These 
two objectives can be met by properly choosing the 
functions A(r) and B(e)« 

Recall that the functions A and B enter 
the transformation through the relationship given 
in Eq. (32) which is now written as 


The derivative E' is 



(38) 


where the minus (plus) sign relates to positive 
(negative) values of r . 


By substituting the expressions for E and 
E' into Eq, (37), we obtain 
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913 =f rl A(r)exp - a o H)l cos (e *) 

\ ' m I / L / 



with the same meaning as above attached to the 
sign notation appearing before the last term. 


As mentioned above, the exponential factor serves 
as a damping function to guarantee that the coor- 
dinate transformations are well behaved; i.e., the 
coordinates should not fold nor should the Jacobian 
of the transformation vanish anywhere within the 
region of interest. In this regard, the quantity 
a 0 is a non-negative constant which increases In 
magnitude as the ratio V/fiR = decreases. For 
a given x, the larger the damping the more the 
e-coordinate is forced to assume a purely 
circumferential direction. 


The periodic significance of z m is now 
determined by introducing Ng as being a posi- 
tive integer and allowing the value of e/r m to 
range from -1/2 Ng to 1/2 Ng; here, Ng signi- 
fies the number of blades in the cascade and may 
be either even or odd. Since the value is 
half the distance between neighboring blades, the 
blades will be located at intervals of twice e m . 
The total distance from -1/2 Ng to 1/2 Ng 
equals one transversal around the cascade. By 
arbitrarily specifying one blade to be positioned 
at e - 0, the blade locations take on the 
following values of e/e m . 

C B (0, ±2, ±4, ±(N b - 1) ; N B odd 

5 m (0, ±2, ±4 ±(Ng - 2), Ng ; N g even 

(40) 


The expression for the metric tensor g 23 
can be written from Appendix A 7 as 


923 


nr r 

or* 




A'(r)BU) A(r)B'(e) 


(45) 


From Eq. (35) and the values of eg* the term 
A*(r)B(c) = 0 at the Ng blade locations. 
Furthermore, since g ^3 is zero at these 
positions, inspection of Eq. (36) indicates that 
AB'*= (nR/V)(r/R) 2 at the blades also. Then, 
for this case, the metric tensor g 23 reduces 
to 


g 23 




fir r IT 
U R U 


/vu ' V'\ nR /r\ 2 

(y--rj Y T w 


where the subscript B denotes a blade location. 
Also, at this point, the meaning of tjie over-bar 
above e is made clear by defining e to be 
the value of z as measured from the nearest 
blade. Thus the extremes of z will be 
With these definitions Eq. (39) for the metric 
tensor can be written 



for values at the blade stations^ However, since 
Cb = 0 from the definition of e, and since 
cos(ir SR/Sm) = 1 and sin(ir *>BlZ m ) * the 
value of the metric tensor g^g for the Ng 
blade stations reduces to 


g 13 


e=e 


B 


fir r V 
= TT R “ H 


A(r) 


(42) 


Therefore, gx 3 can be forced to equal zero at 
these periodic positions by simply prescribing 
that 


fiR 


= ~ V 



(46) 


Requiring V to be a constant is sufficient to 
render g 23 = 0 since then V* * 0. We are in 
this way able to make g £3 as well as gi 3 
equal to zero at the blade stations defined in 
Eq. (40). 

Thus the ^-coordinate lines will be orthogo- 
nal to the helical sheets containing the mean camber 
locations of the blades provided that A(r) and 
B(e) are defined as above and that the axial 
velocity— V — ts_a constant. This is no great 
restriction since a constant value of V is the 
most reasonable case and the one of most interest. 


The second objective is the requirement that the 
axial distance be the same for corresponding points 
at each periodic boundary. This is easily verified 
as being satisfied by inspecting the transformation 
of z on the helical sheets lying midway between 
successive blades. From Eq. (44) it is seen that 


Z PB = 



V 

Y U 


(47) 


when 


fiR 


A(r) = ^ ft 


(43) 


The expression for the axial coordinate can 
now be explicitly written as 

■- -*■ 14,1 


*1, *3 *(N b - 2), N b ; N r odd 

*1, ±3,...,*(N B - 1); Ng even 

where the subscript PB denotes “periodic bound- 
ary." This relation shows that the net axial dis- 
tance is not changed upon complete traversal via 
a e-coordinate line from any periodic sheet to 
any other. In fact, moving a value of twice r m 
along any c-coordinate restores the original 
axial location. That this is true can be seen by 
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inspecting Eq. (44). The consequence of this is 
that any set of sheets separated by a value equal 
to the blade spacing 2 z m could be used as 
periodic boundary sheets.' 

Final Form of the Periodic Helical Coordinates 


For convenience, the transformation of the 
helical coordinates y = y , r = y , I = y 3 to the 
Cartesian coordinates x = x 1 , y = xS z = x 3 is 
given in terms of the nonsuperscripted variables 



This is the set of coordinates which is used in 
the solution of the small disturbance equation. 

Various other cases of these coordinates could 
have been specified. For example, if either A(r) 
or B(0 of Eq. (32) were chosen equal to zero, 
then the t-curves become circles about the 
z-axis. As another special case, if B(r) = z is 
chosen, then the s -curves are helices. On each 
C-helix the coordinate value varies as the arc 
length divided by the value of the radius at which 
that given helix lies. This distinguishes them 
from the streamwise helices which measure as the 
arc length. The metric tensors are given in the 
latter part of Appendix A' for these special 
cases. 


Numerical Approach 


This section presents a general description of 
the numerical approach used to solve the small 
disturbance Eq. (38). For the sake of conven- 
ience, the small disturbance equation will be 
rewritten letting y* = y, y 2 = r, and y 3 = c, as 
in the last section. With this notation Eq. (38) 
can be written 


A. 


A 

3 <p 

1 3 y3 1 


f- L(!*Y + + A- 

3 y 2 W/ 3 3 yJ 4 ar 2 


+ A, 


+ A, 


3 2 CP 


+ A 


+ A 0 


6 3 ysr 7 3 Y 3 C 8 3 r 3 C 


5 77 


(51) 


where 


2MM f c 

A 1 55 277 

6 

A 2 = - \ (/+ ljMMj-g 11 (52b) 

11 m2 

< 52c > 

o 

A 4 = “Z73 g22 = 773 (52d) 

0 o 


*5 - 

(52e) 

2e 12 

6 - 771 = 

(52f) 

CO 

1—1 

il 

r-* 

(52g) 

c 23 

k 8 = ~T73 9 
6 

(52h) 


In transferring the two terms on the left-hand 
side of Eq. (38) to the right-hand side in 
Eq. (51), they have been lumped together with the 
understanding that the variation of g 11 with 
respect to y is small; recall that M is only 
a function of r and, therefore, the only depend- 
ence of A? or A 3 on y is through g 3 *. 


It is convenient to define the two bracketed 
terms above by F so that 


= f- [a? (¥) + A , ri 

3 Y ^ V 3 Y / 3 3y I 


(53) 


Equation (51) is then written as 


2 2 2 2 2 
A 3 2 - F + A -5-2- + A -5-2+fl 3 2 + a — 
A 1 3Y3t - Y ft 4 2 ft 5 2 0 3 Y 3r A 7 3 Y 3f 

3 1 os 


+ a 

O 3r3C 


(54) 


Except for the addition of the last two terms 
above, this equation has the same form as the 
equation solved in Isom 5 where an ADI method 
based on the Douglas-Gunn algorithm was used to 
solve the finite difference form of the equation 
for flow about helicopter rotors. The additional 
cross-derivative terms will be handled by general- 
izing the Douglas-Gunn algorithm. With modest 
changes this allowed the numerical code developed 
for the helicopter problem to be used for the 
present work. 

It should be pointed out that Eq. (54) remains 
expressed in terms of unstretched physical vari- 
ables. No mapping of the coordinate system has 
been carried out so as to produce a nonuniform 
grid. For simplicity, the numerical algorithm 
will be presented for the case of no coordinate 
stretching. Following this, the method of intro- 
ducing coordinate stretching will be explained. 


ADI Douglas-Gunn Algorithm 

For three dimensions the ADI technique involves 
splitting the given equation into three separate 
finite difference equations which can be solved 
successively to complete one time-step increment. 

A current estimate of cp, say, tp n is advanced 
to cp n+ l through two intermediate values, which 
will be denoted by 0 * and cp**, to complete 
a single stage of iteration. 

To begin an iteration at time step n, (P n is 
advanced to cp* by solving the first equation 
with only the Y -direction being differenced 
implicitly; this is called the Y-sweep. Next, in 
the r-sweep, <p* is advanced to «>** by solving 
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the second equation with only the r-direction 
being implicit. Finally, in the *;-sweep, ?** is 
advanced to <p n 1 by solving the third equation 
with the s-direction being the only implicit 
direction. 

The three equations for the respective sweeps 
are given for a uniform grid as 

y-SWEEP: 


F 1+1/2 “ A 2 (Vi,J,k)(v!,J,k) 

+ \ a 3 s y (' p i ,o,k + (58) 

to be the flux at the midpoint of the I™— cel 1 
on any node line given by J and K. The Murman- 
Cole type-dependent differencing scheme is intro- 
duced to provide stable differencing by defining 


II « Y ( f * - <p n ) = ♦ A 4 « rr , n + A 5#e / 

+ A 6 « Yr <p n + A 7 « yC ," + A 8 « rt / (55a) 

r-SWEEP: 


D y F I " 1 7 [ e l( F I+l/2 " F I— 1 / 2 ) + l 1 “ £ I-l) 

x ( F I-l/2 " F I-3/2)] ^ 59) 

where 


Tt « y U** - ,") = D^Fj ♦ ^ « rr (<p** + /) ♦ A 5 « rF , n 




yr 


8 rr 


(55b) 


and 


e I 


v c > 0 
0, V c < 0 


(60) 


e-SWEEP: 

A 1 . , n+1 n. n P . A 4 , ** . rn 

It S y ( <p ~ B °/i 2“ { rr * ) 


+ A 5 , n+1 n> A 6 4 , ** . n, 

T“ ^ + v + (<P + <P ) 


A A 

,7 , n + l . n Y . 8 # , n+1 , n x /cc \ 

y (<p + <p ) + y (<p + <p ) (55c) 


where D y is a special difference operator to be 
explained below and 6 rr , fij-r,..., 6 re are standard, 
difference operators in the indicated directions. 
Here, F represents the bracketed terms of Eq. (51) 
which correspond to the flux in the y-direction, 
and D y Fj gives the finite difference approxima- 
tion to aF/ay at the I T ^-node for each set of 
values of 0 and K. The nodal values I, J and 
K are associated with the y, r and s direc- 
tions, respectively. The three directions, have 
uniform step-sizes Ay, Ar and ae, which may 
be distinct. 


An example of how the difference operators in 
Eq. (55) are defined in terms of difference approx- 
imations is 


5 <p 
rr 


Ur)‘ 


(’I. 


J+1,K 


- 2 cp 


I.J.K 



V c = A 2^" A 3 < 61 > 

This switches the difference equations at each 
grid point according to whether the flow field at 
that point is subsonic, sonic, supersonic, or a 
shock. 

A more convenient set of equations for numer- 
ical computations is obtained from the set (55) 
by subtracting Eq. (55a) from Eq. (55b), and 
Eq. (55b) from Eq. (55c), giving, with some 
rearrangement, the following set: 


y-SWEEP: 
A 


It *,(•* - *"> - °/t ■ Vrr>" * Vt/ 


* V,/ * V,/ * Vr/ 


r-SWEEP: 


It S ‘ T s rr “ T \r, 
A, A 


6 * ^(®** - 

44 Vr ■ 4 \r) (,p * ■ q,n) 


(62a) 


(62b) 


t-SWEEP: 


An analogous expression holds for 6^. In the 
case of the mixed second order operators the fol- 
lowing is used: 


5 


n 


yr 


1 

AyAr 


’I.W.K 


n 

•l.J.K 


n + n 

•l-l,J+l,K 



a 


7 A 
— 0 




A 8 

T 


rz 


n+1 


<p 





(57) 

Corresponding expressions hold for and 6 rj r. 

The nonlinear term contained in D y Fi is 
linearized by averaging at the n and * time 
levels by defining 


The Eq. (62) are in the so-called de’ta form 
where the unknowns on the left-hand sides are 
expressed as differences in the potential. This 
form has superior numerical properties as compared 
to Eq. (55). The solution to the above set of 
equations involve no more than solving tridiagonal 
matrices, except in the case of a shock point where 
a quadridi agonal matrix occurs in the y-$weep. 
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The potential at time level n + 1 can be found 
from the potential at level n by adding the solu- 
tions for the delta differences from all sweeps 
to «> n as 


p P 1 = (p P + (cp* - *") + (q>** - <p*) + (</ 


As is common in ADI methods, the time step is 
varied from iteration to iteration over the course 
of the calculation from some maximum value At max 
to some minimum value At ml * n . This is done to 
improve the convergent rate of the calculation. 

The geometric sequence 

i-1 




i s 1,2,3, . . .N (64) 


is used for N iterations (N = 8) and then 
repeated until the total number of iterations has 
been reached. The total number of iterations is 
determined by a preset value for the maximum num- 
ber of iterations, or by either satisfying a con- 
vergence criteria or exceeding an error bound. 

The range in the time step for the sequence of 
iterations addresses both high and low frequency 
components of the error. In general, the range 
of At min and At max must be determined by trial 
and error and is strongly influenced by the size 
of the computational mesh. 


Grid Stretching 


Up until now the mesh has been considered to 
be a uniform grid. However, it is preferable to 
have the grid points clustered in regions of high 
gradients and sparsely distributed in regions of 
low gradients. Grid stretching is a means of 
accomplishing this. It is used here to distribute 
the physical coordinates near the airfoil with the 
greatest concentrations near the leading and trail- 
ing edges. The grid is smoothly stretched from 
the airfoil surfaces to a coarse grid at the outer 
boundaries of the flow field. The stretching is 
performed in all three coordinate directions. It 
is defined in a general sense as a mapping of the 
physical space yri to a uniform computational 
space yff by 


y = y(y,r), r = r(r) , t = r(e) (65) 

It should be noted that the physical coordinates 
r and K are stretched by their respective com- 
putationaj coordinate. However, y is a function 
of both 7 and r. the added dependence of y 
on r is necessary to accommodate swept wings. 

The coordinate stretching is introduced 
through the familiar chain rule formulas. By 
replacing the various partial derivatives in 
Eq. (51) with those obtained by the chain rule, 
this equation is generalized to stretched grids. 

The use of coordinate stretching complicates, but 
does not change, the basic form of the ADI 
algori thm. 

Discussion and Results 

This discussion focuses on the results obtained 
for flow surrounding a single blade of an eight 


bladed cascade. As explained in the introduction, 
the blades of this cascade are simple bicircular 
arc profiles with a maximum blade thickness of 
five percent. Furthermore, the planform of the 
blades is rectangular with the required spanwise 
twist made about the leading edge. This produces 
a blade of constant chord length. The aspect 
ratio of the blades is defined as the ratio of 
the blade-tip radius, as measured from the axis, 
divided by the chord length. All the results 
here are for an aspect ratio of 4:1. The hub of 
the propeller system is placed at a radius of 
0.375 R, where R is the tip radius; this gives 
an effective aspect ratio of 2.5 for a blade 
length measured from the hub, rather than from 
the axis. 

The following flows were computed: (1) high 

advance ratio x = 100 with a helical free-stream 
Mach number at the blade tip of Mr = 0.8 and 
an angle of attack a * 0.0 ; (2) advance ratio 
x = 1 with Mr = 0.8 and a = 0.0°; (3) o advance 
ratio x = 1 with Mr = 1.1 and a = 0.0°, and 
(4) advance ratio x = 1 with Mr = 1.1 and 
a = 2.0°. The first three cases were compared 
against results from an Euler code. 8 

Before discussing the results of these cases, 
the grids used in the computations for the helical 
small disturbance (HSD) code and the Euler code 
will be discussed. 

Mesh lines for the HSD computations are shown 
(for the case x = 1) in Fig. 6 where, for clarity, 
only every third line is included from the leading 
edge to the trailing edge. A uniform grid is used 
over the blade surface, and a stretched grid is 
used in the following regions: from the leading 

edge to the upstream boundary, from the trailing 
edge to the downstream boundary, from the blade 
tip to the outer radial boundary, and from the 
blade surface to the periodic boundary. The grid 
contained 31 points along the blade in the stream- 
wise direction and an additional 11 grid points 
both before and after the blade, for a total of 
53 streamwise points. In the radial direction, 

22 grid points were used with half of these being 
on the blade. In the circumferential direction, 

22 points were used from the lower to the upper 
periodic boundary. The total grid, therefore, 
was 53x22x22. 

For the Euler computation, the same unstretched 
grid was used over the blade surface. Otherwise, 
the grid was stretched as above but with a differ- 
ent stretching function and, no stretching was 
used in the circumferential direction. The grid 
is shown in Fig. 7 for the case of x = 1. The 
position of the blade is indicated by the narrow 
opening (visible in Fig. 7(b)) of the grid lines 
near the hub. The grid of Fig. 7(c)) illustrates 
how the streamwise grid transitions from the axial 
direction upstream of the blade to a helical direc- 
tion and back again to an axial direction down- 
stream. Since the blade-to-blade direction is 
purely circumferential , this results in a high 
degree of coordinate shearing at axial locations 
near the blade. In addition, the chordwise dis- 
tribution of grid points is not symmetric front to 
back along the blade surface, nor is it symmetric 
from the suction to the pressure side; this asym- 
metry increases with blade thickness and stagger. 
The total grid is 59x21x21 for the case of 
x = 100, and 62x21x21 for the case of x * 1. 
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The above grids are coarse and the computa- 
tional solutions given below will be discussed 
with the understanding that these are preliminary 
results which should be followed by more detailed 
calculations. 

Case 1) x = 100, Mp = 0.8, « = 0° 

This case was chosen to examine what effect 
blade cascading has on the solution. For the 
value of x = 100, the flow is essentially an 
axial flow. Since the blade is symmetric from # 
front to back and from top to bottom and a - 0°, 
the solution should reflect this symmetry if no 
losses occur in the flow field. For the value of 
Mr = 0.8 and the thin 5 percent thick blade, no 
shocks occurred in the flow field. This expected 
symmetry is noticable in the solutions of both 
the HSD and the Euler codes. Figure 8 shows HSO 
Mach contours on the blade surface with the mini- 
mum contour being 0.75 and the maximum 0.9; the 
results are identical for the pressure surface 
and the suction surface. The results reveal the 
expected drop-off in Mach number with increasing 
radius. Similar results are shown in Fig. 9 where 
identical contour values are plotted for the Euler 
code. The contour shapes obtained from the two 
codes are very similar with the only essential 
differences being that the Mach number produced at 
a given blade location is higher for the HSD code. 

Similar symmetry is shown for Mach contours in 
cross section planes given in Figs. 10 and 11 for 
the HSD and Euler codes, respectively. Again, the 
shapes are similar between the two codes, with the 
HSD results showing more flow acceleration through 
the passage. The similarity in shapes indicates 
that qualitatively the HSD solution is being cal- 
culated correctly within the interior region of 
the flow as well as at the blade surfaces. 

The solutions on blade-to-bl ade surfaces are 
given in Figs. 12 and 13 for the HSD and Euler 
codes respectively. In each case Mach contours 
are shown for three different span stations along 
the blade. The minimum contour value is 0.8. The 
values of the maximum contours are as follows: 

1.0 for Fig. 12(a); 0.92 for Fig. 12(b); 0.86 for 
Fig. 12(c); 0.88 for Fig. 13(a); and 0.86 for 
Figs. 13(b) and (c). The results are symmetric 
and support the above finding that the HSD code 
predicts flows that agree with the Euler code 
except in magnitude, at least for subsonic axial 
flow. 

Case 2) x = 1, Mr = 0.8, a = 0° 

This case is presented to isolate the effect 
of blade rotation. The free-stream axial Mach 
number is only 0.5657, although Mr = 0.8. The 
effect of operating at a low advance ratio is seen 
in Figs. 14 and 15 which give Mach number contours 
on the blade surface, again for the HSD and Euler 
computations, respectively. These contours are 
given for a range of 0.6 to 0.8 in each case. 
Other than the expected result that the flow Mach 
number would increase toward the blade tip, the 
Euler contours are asymmetric with the pressure 
surface being nearly the inverted image (left to 
right) as compared to the suction surface. This 
might be the result of blade stagger which would 
give an inverted image for a symmetric blade; for 
the case of isentropic flow, the Mach number on 
the pressure surface, at a given chord location. 


would be the same as that on the suction surface 
if its location was measured from the opposite end 
of the blade. To give an example, the maximum 
Mach number might occur at 60 percent chord, for 
a given span station, on the pressure surface; it 
would then have to occur at 40 percent chord on 
the suction surface. The reason for the observed 
difference in the magnitude between the pressure 
and suction contours for the Euler case is not 
known but, may in part be due to the grid asym- 
metry. While the results of the HSD contours are 
shown to be symmetric, there is no reason that 
the maximum Mach number must be at midchord. 

The blade-to-bl ade contours for this case are 
shown in Figs. 16 and 17. The HSD results are 
given for the following range of Mach contours: 

0.6 to 0.66 for Fig. 16(a); 0.66 to 0.71 for 
Fig. 16(b); and 0.75 to 0.78 for Fig. 16(c). The 
Euler results are given for the same respective 
range of Mach contours. The primary difference 
between the two sets of contours is that the Euler 
contours more closely resemble contours about 
isolated blades. In the case of the HSD contours, 
they tend to shift upstream on the pressure side 
and downstream on the suction side of the blade 
so as to gradually join together at midchannel. 

Case 3) x = 1 , Mr = 1.1, a = 0 

This repeats the previous case except that 
now the free-stream Mach number is increased so 
that it has a value of 0.7778 on the axis and a 
helical free-stream value of 1.1 at the blade tip. 
The Mach contours on the blade surface are given 
in Fig. 18 for the HSD computation. The contours 
are shifted toward the trailing edge on both the 
pressure and suction surfaces, which show identi- 
cal contours. Near the tip and trailing edge a 
very weak shock may exist. In the case of the 
Euler computations, the rearward shift of peak 
Mach number is more pronounced. A weak shock 
probably exists on the suction surface where 
larger gradients than on the pressure side are 
indicated. 

The blade-to-bl ade contours for this case are 
shown in Fig. 20 for the HSD computation and in 
Fig. 21 for the Euler computation. For both sets 
of results, the contour levels range from: a) 

0.82 to 1.0; b) 0.87 to 1.06 and, c) 0.99 to 1.16. 
It is not clear that any shock exists for the HSD 
computation. However, a weak shock is observable 
in Fig. 21(c) of the Euler computation; it origi- 
nates near the trailing edge of the suction sur- 
face and extends outward to a position upstream 
of the neighboring blade. 

Case 4) x = 1, Mr = 1.1, a = 2° 

As a final case, the HSD code was used to 
recalculate the above case except that a spanwise 
uniform angle of attack of 2° was used. The Euler 
code was not used for this case. The Mach contour 
plots are given in Fig. 22 for the blade surfaces 
and in Fig. 23 for the blade-to-bl ade surfaces. 

The effect of imposing an angle of attack to the 
blades resulted in a difference between the pres- 
sure and suction contours in the expected direc- 
tion, i.e., the fluid velocity is now higher on 
the suction side. The bl ade-to-bl ade contours 
reveal that larger regions of the fluid are accel- 
erated on the suction side. Although a weak shock 
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may exist on either surface, no shock is noted to 
extend into the fluid from either surface of the 
blade. 

Concluding Remarks and Future Work 

The above test cases indicate that although 
both the HSD and the Euler codes give reasonable 
solutions, enough differences exist between the 
two sets of solutions to warrant further investi- 
gation, Any future effort should focus on resolv- 
ing the difference between the two codes in the 
values of the Mach contours, particularly for 
case 1. Similarly, the difference between the two 
for the blade-to-blade contours, with x = 1, 
needs to be explained. As noted above, the mesh 
used in these calculations was coarse. In addi- 
tion to the coarseness of the mesh, variations in 
mesh characteristics exist between the two. Some 
of the computational differences between the two 
codes may be resolved if a common mesh is used, 
and a finer mesh should provide more detailed and 
accurate computations. As a means of providing a 
common mesh, a new helical coordinate system is 
currently being developed which conforms to the 
blade shape. This means that the small disturb- 
ance boundary conditions will be able to be 
replaced by the actual blade surface boundary 
conditions. The new coordinate system should 
also be readily adaptable to the Euler code. 

This will allow identical grids to be used for 
both sets of computations. An increased effort 
will then be invested into exploring areas such 
as proper specification of the farfield boundary 
conditions, especially those needed for accurate 
unsteady calculations. 
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FIGURE 1. - RELATIONSHIP BETWEEN INERTIAL COORDINATE 
SYSTEM XYZ AND NON INERTIAL COORDINATE SYSTEM XYZ 
WHICH IS ROTATING ABOUT AXIS PASSING THROUGH ITS 
ORIGIN. 


Z 



FIGURE 2. - NONINERTIAL COORDINATE SYSTEM 
XYZ ROTATES WITH ANGULAR VELOCITY 
ABOUT z-AXIS WHICH IS COL I NEAR WITH Z-AXIS 
OF INERTIAL SYSTEM XYZ. THE XYZ SYSTEM 
TRANSLATES WITH AXIAL VELOCITY V WITH RE- 
SPECT TO THE INERTIAL SYSTEM. 



FIGURE 3. - ARRANGEMENT OF HELICAL CO- 
ORDINATES: THE y 1 -CURVES LIE ALONG 
THE MEAN FLOW DIRECTION: THE BLADE 
PROFILE IS GIVEN BY F y ON THE SUC- 
TION SIDE AND BY F q ON THE PRESSURE 
SIDE; THE FLOW TANGENCY BOUNDARY CON- 
DITION IS NOT APPLIED AT THE BLADE 
SURFACE BUT RATHER ON TWO COORDINATE 
SURFACES NEIGHBORING THE MEAN BLADE 
LOCATION. 



FIGURE A. - PERIODICAL HELICAL COORDINATES 
Yr£ SHOWN WITH CARTESIAN COORDINATES xyz. 
THE y -CURVES ARE ESSENTIALLY RADIALLY DI- 
RECTED; THE CURVES ARE ESSENTIALLY CIRCUM- 
FERENTIAL BUT BEND SO AS TO BE NORMAL TO 
THE y-CURVES AT THE BLADE LOCATIONS. THE 
ORIGIN OF BOTH SYSTEMS IS ON THE AXIS. 



FIGURE 5. - ALTERNATING ARRANGEMENT OF HELI- 
COIDAL SHEETS FOR AN EIGHT BLADED PROPELLER 
PERIODIC BOUNDARY CONDITIONS ARE APPLIED 
AT THE MIDCHANNEL HELICAL SHEETS. THE 
ANGULAR DISTANCE BETWEEN THE TWO TYPES OF 
SHEETS IS DENOTED BY 
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(A) VIEW NORMAL TO AXIS SHOWING PROJECTED 
VIEW OF HELICAL SHEET CONTAINING BLADE. 



(B) CROSS-CHANNEL VIEW AT LEADING EDGE. 



(C) BLADE-TO-BLADE VIEW 
ACROSS TWO CHANNELS; 

NOT ALL GRID POINTS ARE 
SHOWN ALONG BLADE SURFACE. 


FIGURE 6. - HELICAL COORDINATES USED IN SMALL DISTURBANCE COMPUTATION. 




(B) VIEW ALONG AXIS SHOWING THE MIDCHORD 
STATION. 


(C) VIEW OF BLADE-TO-BLADE SURFACE 
SHOWING THREE BLADES. GRID SHOWN 
IS COARSER IN AXIAL DIRECTION THAN 
THAT USED. 


FIGURE 7. - COORDINATES USED IN EULER COMPUTATION. 
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(A) PRESSURE SURFACE. (B) SUCTION SURFACE. 

FIGURE 8. - MACH CONTOURS OF SMALL DISTURBANCE 
COMPUTATION ON BLADE SURFACES: ADVANCE RATIO = 
RATIO = 100. M r = 0.8, a = 0°. 


(A) PRESSURE SURFACE. <B) SUCTION SURFACE. 

FIGURE 9. - MACH CONTOURS OF EULER COMPUTATION 
ON BLADES SURFACES: ADVANCE RATIO = 100, 

M r = 0.8, a = 0°. 
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HUB TIP 


FIGURE 10. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION IN 
CROSS PLANE AT MIDCHORD AXIAL LOCATION: ADVANCE RATIO = 100, 
M r = 0.8,0 = 0°. 



HUB TIP 


FIGURE 11. - MACH CONTOURS OF EULER COMPUTATION IN CROSS PLANE 
AT MIDCHORD AXIAL LOCATION: ADVANCE RATIO = 100, M D = 0 8, 

a = 0°. R 



(A) r 0.375 R. <B) r - 0.625 R. (C) r = 0 .875 R. 

R m UR - n 2 « ' MACH 0 CONTOURS 0,1 DISTURBANCE COMPUTATION ON BLADE-TO -BLADE SURFACE AT VARIOUS SPAN LOCATIONS: ADVANCE RATIO = 100 

Mb 0.8, G - 0 . 
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(A) PRESSURE SURFACE. 


(B) SUCTION SURFACE. (A) PRESSURE SURFACE. 


(B) SUCTION SURFACE. 


FIGURE 14. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION ON FIGURE 15. - MACH CONTOURS OF EULER COMPUTATION ON BLADE SURFACES: 

BLADE SURFACES: ADVANCE RATIO = 1, M R = 0.8, a = 0°. ADVANCE RATIO = 1, M R = 0.8, a = 0°. 
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(C) r = 0.875 R. 


FIGURE 16. - MACH CONTOURS OF SMALL DISTURBANCE COMPUTATION ON BLADE -TO-BLADE SURFACE AT VARIOUS SPAN LOCATIONS; ADVANCE RATIO = 1, 
M r = 0.8, a = 0°. 
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(B) r = 0.625 R. 
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(C) r * 0.875 R. 


FIGURE 17. - MACH CONTOURS OF EULER COMPUTATION ON BLADE-TO -BLADE SURFACE AT VARIOUS SPAN LOCATIONS; ADVANCE RATIO - 1, M R = 0.8, 
a = 0°. 
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(A) r = 0.375 R. 
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(B) r = 0.625 R. 
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(C) r = 0.875 R. 


FIGURE 21. - MACH CONTOURS OF EULER COMPUTATION ON BLADE-TO-BLADE SURFACE AT VARIOUS SPAN LOCATIONS; ADVANCE RATIO = 1. It * 1.1. 
a = oo_ 
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(A) PRESSURE SURFACE. (B) SUCTION SURFACE. 

FIGURE 22. - MACH CONTOURS OF SMALL DISTURBANCE ON BLADE SUR- 
FACES: ADVANCE RATIO = 1 , M R = 1 . 1 , Q s 2° . 
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